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Abstract 

We discuss 3D global simulations of the early Martian climate that we have per- 
formed assuming a faint young Sun and denser CO2 atmosphere. We include a 
self-consistent representation of the water cycle, with atmosphere-surface interac- 
tions, atmospheric transport, and the radiative effects of CO2 and H2O gas and 
clouds taken into account. We find that for atmospheric pressures greater than a 
fraction of a bar, the adiabatic cooling effect causes temperatures in the south- 
ern highland valley network regions to fall significantly below the global average. 
Long-term climate evolution simulations indicate that in these circumstances, wa- 
ter ice is transported to the highlands from low-lying regions for a wide range of 
orbital obliquities, regardless of the extent of the Tharsis bulge. In addition, an 
extended water ice cap forms on the southern pole, approximately corresponding to 
the location of the Noachian/Hesperian era Dorsa Argentea Formation. Even for a 
multiple-bar CO2 atmosphere, conditions are too cold to allow long-term surface liq- 
uid water. Limited melting occurs on warm summer days in some locations, but only 
for surface albedo and thermal inertia conditions that may be unrealistic for water 
ice. Nonetheless, meteorite impacts and volcanism could potentially cause intense 
episodic melting under such conditions. Because ice migration to higher altitudes 
is a robust mechanism for recharging highland water sources after such events, we 
suggest that this globally sub-zero, 'icy highlands' scenario for the late Noachian 
climate may be sufficient to explain most of the fluvial geology without the need to 
invoke additional long-term warming mechanisms or an early warm, wet Mars. 
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1 Introduction 



After many decades of observational and theoretical research, the nature of 
the early Martian climate remains an essentially unsolved problem. Extensive 
geological evidence indicates that there was both flowing liquid water (e.g., 



Carr 


|1995 , 


Irwin et al. 


2005], 1 



and standing bodies of water (e.g., Fassett and Head [2008bj) on the Martian 



surface in the late Noachian, but a comprehensive, integrated explanation for 
the observations remains elusive. As the young Sun was fainter by around 25 
% in the Noachian (before approx. 3.5 GYa) [Gough , 1981 , to date no climate 



model has been able to produce long-term warm, wet conditions in this period 
convincingly. Transient warming events have been proposed to explain some of 
the observations, but there is still no consensus as to their rate of occurrence 
or overall importance. 



The geomorphological evidence for an altered climate on early Mars includes 
extensive dendritic channels across the highland Noachian terrain (the famous 
Valley networks') |Carr[ |1996[ [Fassett and Headj |2008b[ [Hynek et"aL| |2010 



et al. 



fossilized river deltas with meandering features [Malin and Edgett[ [2003 , Fas 
2005|, records of quasi-periodic sediment deposition [Lewis 



sett and Head 

2008], and s of enhanced erosion most readily explained through 



fluvial activity [Hynek and Phillips , 2001J . Some studies have also suggested 
evidence for an ancient ocean in the low- lying northern plains. These include 
a global analysis of the martian hydrosphere [Clifford and Parker, 2001 and 
an assessment of river delta / valley network contact altitudes |di Achille and 
Hynek, 2010 . However, in the absence of other evidence, the existence of a 



northern ocean in the Noachian remains highly controversial. 

More recent geochemical evidence of aqueous alteration on Mars has both 
broadened and complicated our view of the early climate. Observations by the 
OMEGA and CRISM instruments on the Mars Express / Mars Reconnais- 
sance Orbiter spacecraft iPoulet et al., 2005, Bibring et al., 2006, Mustard 



et al., 2008, Ehlmann et al., 2011 showed widespread evidence for phyllosili- 
cate (^clay) and sulphate minerals across the central and southern Noachian 
terrain. Surface aqueous minerals are rarer in Mars' northern lowlands, which 



are mostly covered by younger Hesperian-era lava plain s ^Head et al. , 2002 



Salvatore et al., 2010[ and outflow channel eflluent [Kreslavsky and Head 
2002[. However, phyllosilicates have been detected in some large northern im 



pact craters that penetrated through these later deposits Carter et al., 2010 
As these impacts were understood to have excavated ancient Noachian terrain 
from below the lava plains, it seems plausible that aqueous alteration was once 



^ Current address: Department of the Geological Sciences, University of Chicago, 
60637 IL, USA. 
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widespread in both hemispheres, on or just beneath the Martian surface. 



Evidence from crater counting |Fassett and Head , 2008a, 2011| suggests that 
valley network formation was active during the Noachian but ended quite 
sharply at the Noachian-Hesperian boundary (approx. 3.5 GYa in their anal- 
ysis). Broadly speaking, this period overlaps with the end of the period when 
impacts were frequent and the Tharsis rise was still forming. Interestingly, 
however, crater statistics also suggest that the main period of phyllosilicate 
formation ended somewhat before the last valley networks were created. Few 



Late Noachian open-basin lakes Fassett and Head, 2008b show evidence of 
extensive in-situ phyllosilicates on their floors Goudge et al., 2012 and in 



those that do, the clays appear to have been transported there from older 



deposits by way of valley networks [e.g., Ehlmann et al., 2008j 



Interpreting the surface conditions necessary to form the observed phyllosil- 
icates on Mars remains a key challenge to understanding the Noachian cli- 
mate. It is clear that there is substantial diversity in the early Martian miner- 
alogical record, which probably at least partially reflects progressive changes 



in environmental conditions over time [Bibring et al.. 


2006 


, Mustard et al.. 


2008, 


Murchie et al., 2009, Andrews-Hanna et al. 


, 2010 


Andrews-Hanna and 


Lewis 


, 2011 . Nonetheless, the most recent reviews of t' 


le available geochem- 



ical evidence suggest that the majority of phyllosilicate deposits may have 
been formed via subsurface hydrothermal alteration | Ehlmann et al., 2011 or 



episodic processes, as opposed to a long-term, warm wet climate. 



While it is likely that Mars once possessed a thicker CO2 atmosphere than it 
has today, it is well known that CO2 gaseous absorption alone cannot produce 
a greenhouse eflFect strong enough to allow liquid water on early Mars at any 



atmospheric pressure [Kasting , 1991 . Various alternative explanations for an 
early warm, wet climate have been put forward. Two of the most notable are 
additional absorption by volcanically emitted sulphur dioxide [Halevy et al. , 
20071, and downward scattering of outgoing infrared radiation by CO2 clouds 



Forget and Pierrehumbert , 1997|. However, both these hypotheses have been 
criticised as insufl&cient in later studies [Colaprete and To"on , 2003, Tian et al. 



2010| . Sulphur-induced warming is attractive due to the correlation between 
the timing of Tharsis formation and the valley networks and the abundance 



2010 



ar- 



of sulphate minerals on the Martian surface. However, Tian et al. 
gued that this mechanism would be ineffective on timescales longer than a few 
months due to the cooling effects of sulphate aerosol formation in the high at- 
mosphere. CO2 clouds are a robust feature of cold CO2 atmospheres that have 



already been observed in the mesosphere of present-day Mars Montmessin 



et al., 2007| . They can cause extremely effective warming via infrared scat- 
tering if they form at the correct altitude and have global coverage close to 



100 %. However, our 3D simulations of dry CO2 atmospheres Forget et al. 



2012| suggest that their warming effect is unlikely to be strong enough to raise 
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global mean temperatures above the melting point of water for reasonable at- 
mospheric pressures. 



Given the problems with steady-state warm, wet models, other researchers 
have proposed that extreme events such as meteorite impacts could be capa- 



ble of causing enough warming to explain the observed erosion alone Segura 



et al. , 


2002, 


2008, 


Toon et al. , 



Toon et al., 2010| . These authors proposed that transient 



steam atmospheres would form for up to several millenia as a result of im- 
pacts between 30 and 250 km in diameter. They argued that the enhanced 
precipitation rates under such conditions would be sufficient to carve valley 
networks similar to those observed on Mars, and hence that a long-term warm 
climate was not necessary to explain the geological evidence. This hypothesis 
has been questioned by later studies - for example, landform evolution mod- 
elling of the Parana Valles region (-20° N, 15° W) [Barnhart et"aL| |2009| has 



suggested that the near-absence of crater rim breaches there is indicative of a 
long-term, semi-arid climate, as opposed to intermittent catastrophic deluges. 
In addition, other authors have noted that there is a significant discrepancy 
(as large as 10"^) between the estimated Noachian erosion rates and the to- 



tal erosion possible due to post-impact rainfall [Ramirez et al., 2012 . Hence 



impact-generated steam atmospheres alone still appear unable to explain key 
elements of the geological observations, and the role of impacts in the Noachian 
hydrological cycle in general remains unclear. 

Most previous theoretical studies of the early Martian climate have used one- 
dimensional, globally averaged models. While such models have the advantage 
of allowing a simple and rapid assessment of warming for a given atmosphere, 
they are incapable of addressing the infiuence of seasonal and topographic 



temperature variations on the global water cycle. Johnson et al. |2008| ex- 
amined the impact of sulphur volatiles on climate in a 3D general circulation 
model (GCM), but they did not include a dynamic water cycle or the radia- 
tive effects of clouds or aerosols. To our knowledge, no other study has yet 
attempted to model the primitive Martian climate in a 3D GCM. 

Here we describe a range of three-dimensional simulations we have performed 
to investigate possible climate scenarios on early Mars. Our approach has 
been to study only the simplest possible atmospheric compositions, but to 
treat all physical processes as accurately as possible. We modelled the early 
Martian climate in 3D under a denser CO2 atmosphere with a) dynamical 
representation of cloud formation and radiative effects (CO2 and H2O) b) 
self-consistent, integrated representation of the water cycle and c) accurate 
parameterisation of dense CO2 radiative transfer. We have studied the effects 
of varying atmospheric pressure, orbital obliquity, surface topography and 
starting H2O inventory. In a companion paper fForget et all 120121 , we describe 



the climate under dry (pure CO2) conditions. Here we focus on the water cycle, 
including its effects on global climate and long-term stabilization. Based on 
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our results, we propose a new hypothesis for vaUey network formation that 
combines aspects of previous steady-state and transient warming theories. 



In Section 2, we describe our chmate model, including the radiative transfer 
and dynamical modules and assumptions on the water cycle and cloud forma- 
tion. In Section 3, we describe the results. First, 100% relative humidity sim- 
ulations are analysed and compared with results assuming a dry atmosphere 
[Forget et al. , 2012 . Next, simulations with a self-consistent water cycle and 
varying assumptions on the initial CO2 and H2O inventories and surface to- 
pography are described. Particular emphasis is placed on a) the long-term 
evolution of the global hydrology towards a steady state and b) local melting 
due to short-term transient heating events. In Section 4 we discuss our results 
in the context of constraints from geological observations and atmospheric 
evolution theory, and assess the probable effects of impacts during a period of 
higher flux. Finally, we describe what we view as the most likely scenario for 
valley network formation in the late Noachian and suggest a few directions for 
future study. 



2 Method 



To produce our results we used the LMD Generic Climate Model, a new cli- 
mate simulator with generalised radiative transfer and cloud physics that we 



Selsis et al. 




2011 


, Wordsworth et al., 


2011 . The model uses 


dynamical core fHourdin et al., 


2006| 


which is based on a 



formulation of the classical primitive equations of meteorology. The numerical 
scheme is formulated to conserve both potential enstrophy and total angular 
momentum [Sadourny , 1975| . 



Scale-selective hyperdiffusion was used in the horizontal plane for stability, 
and linear damping was applied in the topmost levels to eliminate the artiflcial 
reflection of vertically propagating waves. The planetary boundary layer was 



parameterised using the method of Mellor and Yamada [1982 and Galperin 
et al. [1988 to calculate turbulent mixing, with the latent heat of H2O also 
taken into account in the surface temperature calculations when surface ice 
/ water was present. As in Wordsworth et al. [2011 , a standard roughness 



coeflicient of = 1 x 10~^ m was used for both rocky and ice / water surfaces. 
A spatial resolution of 32 x 32 x 15 in longitude, latitude and altitude was used 
for the simulations. This was slightly greater than the 32 x 24 x 15 used in most 
simulations in Forget et al. |2012) to allow more accurate representation of the 
latitudinal transport of water vapour, but still low enough to allow long-term 
simulation of system evolution in reasonable computation times. 
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2.1 Radiative transfer 



Our radiative transfer scheme was similar to that used in previous studies 
[e.g., Wordsworth et al., 2011). For a given mixture of atmospheric gases, 



we computed high resolution spectra over a range of temperatures, pres- 
sures and gas mixing ratios using the HITRAN 2008 database Rothman 
et al., [2009 1 . For this study we used a 6 x 9 x 7 temperature, pressure 
and H2O volume mixing ratio grid with values T = {100, 150, . . . ,350} K, 
p = {10-^10-^...,10^} mbar and qh^o = {10"^ 10"^ . . . , lO"^}, respec- 
tively. The correlated-/^ method was used to produce a smaller database of 
coefficients suitable for fast calculation in a GCM. The model used 32 spec- 
tral bands in the longwave and 36 in the shortwave, and sixteen points for 
the g'-space integration, where g is the cumulated distribution function of the 
absorption data for each band. CO2 collision-induced absorption was included 
using a parameterisation based on the most recent theoretical and experimen- 
tal studies Wordsworth et al., 2010a, Gruszka and Borysow, 1998, Baranov 



et al., 2004 



H2O lines were truncated at 25 cm~ , while the H2O contin- 
uum was included using the CKD model [Clough et al. , p!989 j . In clear-sky 
one-dimensional radiative-convective tests (results not shown) the H2O con- 
tinuum was found to cause an increase of less than 1 K in the global mean 



surface temperatures under early Mars conditions. A two-stream scheme Toon 



et al. , 1989J was used to account for the radiative effects of aerosols (here, CO2 



and H2O clouds) and Rayleigh scattering, as in [Wordsworth et al.| [2010^ 



2.2 Water cycle and CO2 clouds 



Three tracer species were used in the simulations: CO2 ice, H2O ice and H2O 
vapour. Tracers were advected in the atmosphere, mixed by turbulence and 
convection, and subjected to changes due to sublimation / evaporation and 
condensation and interaction with the surface. For both gases, condensation 
was assumed to occur when the atmospheric temperature dropped below the 
saturation temperature. Local mean CO2 and H2O cloud particle sizes were 
determined from the amount of condensed material and the number density 
of cloud condensation nuclei [CCN]. The latter parameter was taken to be 
constant everywhere in the atmosphere. Its effect on cloud radiative effects 
and hence on climate is discussed in detail in Forget et al. [2012 |; here it was 



taken to be a global constant (10 kg ; see Table 1). 

Ice particles of both species were sedimented according to a version of Stokes' 
law appropriate for a wide range of pressures f Rossow , 1978 . Below the strato- 



sphere, adjustment was used to relax temperatures due to convection and / 
or condensation of CO2 and H2O. For H2O, moist and large-scale convection 
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were taken into account following Manabe and Wetherald 1967| . Precipita- 



tion of H2O was parameterized using a simple cloud water content threshold 



Iq Emanuel and Ivkovi-Rothman, 1999 , with precipitation evaporation also 



taken into account. To test the sensitivity to precipitation assumptions, we 
performed some simulations a) using a more complex scheme originally de- 
signed for use on the Earth |Le Trent and Li, 1991| and b) with precipitation 



neglected entirely. In the former case, we found overall similar behaviour, with 
slight variations in the mean water cloud column amount causing a 2 K differ- 
ence in the surface temperature at 1 bar surface pressure. When precipitation 
was removed entirely, we found that much thicker H2O ice clouds formed, 
leading to a surface warming of over 10 K. As particle coagulation is a fun- 
damental physical process in water liquid/ice clouds, we did not regard these 
latter tests as physically robust, and used the simplified threshold scheme for 
all simulations with a water cycle presented here. 

On the surface, the local albedo was varied according to the composition 
(rocky, ocean, CO2 or H2O ice; see Table 1). When the surface ice amount 
Qice was lower than a threshold value Qcovered = 33 kg m~^ (corresponding to 
a uniform layer of approx. 3.5 cm thickness; Le Trent and Li| |1991| ), partial 



coverage was assumed and the albedo Ag was calculated from the rock and 
water ice albedos Aj. and Ai as 



M9,\)=Ar + iAi-Ar)^^ (1) 

Qcovered 



where 9 is longitude and A is latitude. For Qice > Qcovered^ the local albedo was 
set to Ai. When CO2 ice formed, the local albedo was set to Ai immediately. 
Water ice formation (melting) was assumed to occur when the surface tem- 
perature was lower (higher) than 273 K, and the effects of the latent heat of 
fusion on temperature were taken into account. As in most simulations only 
transient melting occurred in localised regions, the effects of runoff were not 
taken into account. Surface temperatures were computed from radiative, sen- 
sible and latent surface heat fiuxes using an 18-level model of thermal diffusion 



in the soil Forget et al., 1999 and assuming a homogeneous thermal inertia 
(250 J m~^ s~^^'^ K"-*^). The dependence of the results on the assumption of 
constant thermal inertia are discussed in Section 3. 



2.3 Ice equilibration algorithm 



The time taken for an atmosphere to reach thermal equilibrium can be esti- 
mated from the radiative relaxation timescale [Goody and Yung" , 1989] 
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where Cp^ g, Te and a are the specific heat capacity of the atmosphere, the 
mean surface pressure, the surface gravity, the atmospheric emission temper- 
ature and the Stefan-Boltzmann constant, respectively. Taking Tg = 200 K, 
Cp = 800 J kg"^, = 1 bar and g = 3.72 m s~^ yields ^ 1 Mars year. 
For simulations without a water cycle, this was a good approximation to the 
equilibration time for the entire climate. 

When a water cycle is present, equilibration times can be far longer. Particu- 
larly for cold climates, where sublimation and light snowfall are the dominant 
forms of water transport, the surface H2O distribution can take thousands of 
years or more to reach a steady state. Running the 3D model took around 
5 hours for one simulated Mars year at the chosen resolution, so evaluating 
climate evolution in the standard mode was prohibitively time-consuming. 

To resolve this problem, we implemented an ice iteration scheme. Starting from 
the initial surface ice distribution, we ran the GCM normally in intervals of 
two years. In the second year, we evaluated the annual mean ice rate of change 
{dhice/dt) at each gridpoint. This rate was then multiplied by a multiple-year 
timestep At to give the updated surface ice distribution 



In addition, all ice was removed in regions where h^^^ dropped below zero or 
ice coverage was seasonal (i.e. hice = at some point during the year). After 
redistribution, the amount in each cell was normalised to conserve the total 
ice mass in the system Mh20- The GCM was then run again for two years and 
the process repeated until a steady state was achieved. 

Trial and error showed that correct choice of the timestep At was important: 
when it was too high, the surface ice distribution tended to fluctuate and did 
not reach a steady state, while when it was too low the simulations took a 
prohibitively long time to converge. We used a variable timestep to produce 
the results described here. For the first 5 iterations. At was set to 100 years, 
after which it was reduced to 10 years for the final 15 iterations. This allowed 
us to access the final state of the climate system after reasonable computation 
times without compromising the accuracy of the results. 

In all simulations we set the total ice mass Mh20 = 4:^R^Qcovered^ with R the 
planetary radius. This quantity was chosen so that for completely even ice 
coverage, the entire surface would have the maximum albedo Ai [see equation 
(1)]. While the total Martian H2O inventory in the Noachian is likely to have 
been significantly greater than this, such an approach allowed us to study the 
infiuence of topography and climate on the steady state of the system without 
using unreasonably long iteration times. Two types of initial conditions were 




(3) 



8 



used: in the first, ice was restricted to altitudes lower than -4 km from the 
datum ('icy lowlands'), while in the second, ice was restricted to latitudes of 
magnitude greater than 80° ('icy poles'). Varying the initial conditions in this 
way allowed us to study the uniqueness of climate equilibria reached using 
equation (3). 



2.J^ Topography 



Techniques such as spherical elastic shell modelling and statistical crater anal- 
ysis date the formation of the majority of the Tharsis rise to the mid to late 
Noachian [Philhps et al. , 2001, Fassett and Head, 2011 . Here, we performed 
most simulations with present-day surface topography, but we also investi- 
gated the climatic effects of a reduced Tharsis bulge. For these cases, we used 
the formula 



^mod{0, A) = - (0 + 4000^)exp[(x/65°) 



o\4.5l 



(4) 



to convert the present-day geopotential to (pmod^ where X = ([^ + 100°]^ + 
A^)°-^ and g is surface gravity. Figure 1 compares contour plots of the standard 
topography and that described by equation (4). 

Table 1 

Standard parameters used in the climate simulations. 



Parameter 


Values 


Orbital eccentricity e 


0.0 


Obliquity [degrees] 


25.0, 45.0 


Surface albedo (rocky) Ar 


0.2 


Surface albedo (liquid water) 


0.07 


Surface albedo (CO2 / H2O ice) Ai 


0.5 


Surface topography 


present-day, no Tharsis 


Initial CO2 partial pressure PCO2 [bars] 


0.008 to 2 


Surface roughness coefficient zq [m] 


1 X 10-2 


Surf, therm, inertia I [J m~^ s~-^/^ K~^] 


250 


H2O precipitation threshold Iq [kg kg~^] 


0.001 


No. of cloud condens. nuclei [CCN] [kg~^] 


1 X 10^ 
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3 Results 



3.1 Fixed atmospheric humidity simulations 



In [Forget et~aL] |2012 , we described the climate of early Mars with a dense, 



dry (pure CO2) atmosphere. Here we begin by considering the effects of water 
vapour on surface temperatures. We performed simulations as in the baseline 
cases of [Forget eta'Ll|2012] (variable CO2 pressure, circular orbit, 25° obliquity. 



present-day topography) for a) pure CO2 atmospheres and b) atmospheres 
with relative humidity of 1.0 at all altitudes. In these initial simulations, wa- 
ter vapour was only included in the atmospheric radiative transfer, and not 
treated as a dynamical tracer or used in convective lapse rate calculations. 
Neglecting water vapour in lapse rate calculations should only create a small 
positive surface temperature error in cold climates, and the effect of water 
clouds on climate in our simulations with a full water cycle was small (this 
is discussed further later). This approach hence allowed us to place an upper 
bound on the possible greenhouse warming in our calculations, independent 
of uncertainties in H2O convection par ameterisat ions. 

In Figure 2, global mean and annual maximum/minimum temperatures from 
simulations with surface pressures ranging from 0.2 to 2 bar are plotted. Com- 
parison with results for clear CO2 atmospheres shows the warming effect of 
CO2 clouds (see Fig. 1 of Forget et al. |2012| ). As expected, in atmospheres 



saturated with water vapour the net warming is greater (by a factor of a few 
Kelvin to ^20 K, depending on the surface pressure). Although there are no 
scenarios where the mean surface temperature exceeds 273 K, annual maxi- 
mum temperatures are significantly greater than this in all cases. For pressures 
greater than 2 bar or less than 0.5 bar, permanent CO2 ice caps appeared in 
the simulations, indicating the onset of atmospheric collapse. As discussed in 
Forget et al. [2012 , finding the equilibrium state in such cases would require 



extremely long iteration times. Between 0.5 and 2 bar, however, CO2 surface 
ice was seasonal, and the climate hence reached a steady state on a timescale 
of order (in the absence of a full water cycle). 

In Figure 3, annual mean (top row), diurnally averaged annual maximum (sea- 
sonal maximum; middle row) and absolute annual maximum (bottom row) 
surface temperatures in 2D are plotted for atmospheres of 0.008 (left), 0.2 
(middle) and 1 (right) bar pressure in the maximally H20-saturated simu- 
lations in equilibrium. The diurnal averaging consisted of a running one-day 
mean over results sampled 8 times per day. The adiabatic cooling effect de- 
scribed in Forget et al. |2012| is clearly apparent in the annual mean temper- 



atures (Figure 3, top row): at 0.008 and 0.2 bar (left and middle) the main 
temperature difference is between the poles and the equator, while at 1 bar 
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(right) altitude-temperature correlations dominate. With an atmospheric pres- 
sure close to that of Earth, the regions of Mars with the most concentrated 
evidence for flowing liquid water (the southern Noachian highlands) are among 
the coldest on the planet. 

The seasonal and absolute maximum temperatures are partially correlated 
with altitude even at 0.2 bar (Figure 3, middle column, middle and bottom), 
with the result that only small regions of the southern highlands are ever 
warmer than 273 K at that pressure. At 1 bar, seasonal maximum tempera- 
tures (Figure 3, right column, middle) are only above 273 K in the northern 
plains, parts of Arabia Terra and the Hellas and Argyre basins. However, 
the absolute maximum temperatures exceed 273 K across most of the planet 
(Figure 3, right column, bottom) except the Tharsis rise and southern pole. 
This has interesting implications for the transient melting of water ice, as we 
describe in the next section. 



3.2 Simulations with a water cycle 



When the water cycle is treated self-consistently, surface temperatures can dif- 
fer from those described in the last section due to a) variations in atmospheric 
water vapour content, b) the radiative effect of water clouds and c) albedo 
changes due to surface ice and water. We studied the evolution of the surface 
ice distribution under these conditions using the ice evolution algorithm de- 
scribed in Section 2. We tested the effects of obliquity, the surface topography 
and the starting ice distribution in the model. 



In general terms, ice should tend to accumulate over time in regions where it 
is most stable; i.e., the coldest parts of the planet. On present-day Mars, the 
atmospheric influence on the surface energy balance is small, and the annual 



mean surface temperature is primarily a function of latitude [e.g.. Forget et al. 



1999| . However, as discussed in the previous section, increased thermal cou- 



pling between the atmosphere and the surface at higher pressures in the past 
[Forget et al. , 2012 would have caused a signiflcantly greater correlation be- 
tween temperature and altitude, suggesting that ice might have accumulated 
preferentially in highland regions then. In reality, atmospheric dynamics can 
also play an important role on the H2O surface distribution. This can occur 
through processes such as topographic forcing of precipitation or the creation 
of global bands of convergence / divergence (c.f. the inter-tropical convergence 
zone on Earth, e.g. Pierrehumbert 2011| ). The complexity of these effects is an 
important reason why 3D circulation modelling is required for a self-consistent 
analysis of the primitive Martian water cycle. 
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3.2.1 Surface ice evolution 



Figures 4-7 show the evolution of the surface H2O ice distribution for several 
simulations with different initial conditions and topography. Snapshots of the 
ice are given at the start of the 1"*, 4*^^, 20**^ and 40*^^ years of simulation, with 
the iteration algorithm applied as described in Section 2.3. Simulations were 
run for a range of pressures; here we focus on results at 0.008, 0.04 and 1 bar 
only. 

Figure 4 shows two simulations with identical climate parameters (1 bar pres- 
sure, (/) = 25°) but different initial surface distributions of water ice ('icy 
lowlands' and 'icy poles' for left and right columns, respectively). As can be 
seen, the two cases evolve towards similar equilibrium states, with ice present 
at both poles and in the highest-altitude regions across the planet (the Tharsis 
bulge, Olympus and Elysium Mons and the Noachian terrain around Hellas 
basin). The presence of significant amounts of ice over the equatorial regions 
means an increased planetary albedo. As a result of this and the reduced 
relative humidity due to localisation of surface H2O sources, the mean sur- 
face temperature in the final year (^230 K) is around 5 K below that of the 
equivalent 100% humidity simulation. Figure 5 shows the longitudinally and 
yearly averaged surface ice as a function of time and latitude for the same 
two simulations. The plots indicate that even at the end of the simulation, the 
ice density around the equator and in the south was still increasing slightly, 
mainly at the expense of deposits around 40 °S. Analysis of the atmospheric 
dynamics (see next section) suggested that low relative humidity associated 
with convergence of the annual mean meridional circulation at mid-latitudes 
was responsible for this long-term effect. 

In Figure 6 (left column), the same results are plotted for a 1 bar, cj) = 25° case 
with modified surface topography as described by equation (4) . Here evolution 
is similar to that in the standard cases, except that the correlation between ice 



deposits and the distribution of Noachian-era valley networks [see e.g. Fassett 



and Head, 2008b, Hynek et al. , 2010 is even clearer because the Tharsis bulge 
is absent. When the obliquity is increased, results are also broadly similar 
(see (/) = 45° case; right column of Figure 6), except that ice disappears from 
the northern pole entirely due to the increased insolation there. The effects 
of obliquity on ice transport has been discussed in the context of the present- 



day Martian atmosphere by Forget et al. (2006| and Madeleine et al. |2009j . 



In these dense-atmosphere simulations, the combination of adiabatic heating 
and increased insolation in the northern plains at high obliquities makes them 
the warmest regions of the planet, with mean temperatures of ^ 235 K even 
for latitudes north of 5 °N. 

As a test of the ice evolution algorithm, we also performed simulations at 
0.008 and 0.04 bar pressure (Figure 7). Converging on a solution proved chal- 
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lenging in these the rate of H2O ice transport was extremely slow, 

and permanent CO2 ice caps formed at the poles as on present-day Mars. 
At 0.008 bar, the total atmospheric mass was so low that after several years 
the simulations approached a seasonally varying CO2 vapour-pressure equilib- 
rium. At 0.04 bar, however, atmospheric pressure continued to decrease slowly 
throughout the simulation. 

As a result of the large H2O ice evolution timescales, the simulations did not 
reach an exact equilibrium state even after the full 40 years of simulation time. 
Nonetheless, Figure 7 clearly shows the differences from the 1-bar case: H2O 
ice is present in large quantities at the poles, with a smaller deposit over the 
Tharsis bulge. Analysis of the long-term ice trends in these cases (not shown) 
indicated that the southern polar caps were still slowly growing even at the 
end of the simulations. 

The differences in ice migration between the low and high pressure simulation 
are easily understood through the adiabatic effect described earlier and shown 
in Figure 3. In all cases, ice tends to migrate towards the coldest regions on 
the surface. At 0.008 bar (Figure 3, left-hand column), these regions correlate 
almost exactly with latitude, while at 1 bar (Figure 3, right-hand column), 
the correlation is primarily altitude-dependent. 



3.2.2 Atmospheric dynamics and clouds 

Although the primary factor in the H2O ice distribution is the annual mean 
surface temperature distribution, the final steady state was modulated by the 
effects of the global circulation. Figures 8 and 9 show the annual mean hor- 
izontal velocity at 500 hPa (left) and the annual and zonal mean vertical 
velocity cu (right) for the final year of the 1-bar simulation with obliquity 25° 
and 45°, respectively. As can be seen, the annual mean circulation consists of 
an equatorial westward jet, with a transition to eastward jets above absolute 
latitudes of 40°. The associated vertical velocities are asymmetric as a result 
of the planet's north/south topographic dichotomy, with intense vertical shear 
near 20°N in both cases. Nonetheless, wide regions of downwelling (positive cu) 
are apparent around 50°N/S. In analogy with the latitudinal relative humid- 
ity variations due to the meridional overturning circulation on Earth, these 
features explain the tendency of ice to migrate away from mid-latitudes, as 
apparent from Figures 4-6 (described in the last section). 

Figure 10 shows the total precipitation in each season for the same simulation 
(obliquity 25°). It highlights the dramatic annual variations in the meteorol- 
ogy. From Lg = 0° to 180° (top two maps), relatively intense precipitation 
(snowfall) occurred in the northernmost regions due to the increased surface 
temperatures there. In contrast, from = 180° to 360° (bottom two maps) 
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lighter snowfall occurred in the south, particularly in the high Noachian terrain 
around Hellas basin where the valley network density is greatest. 

Figure 11 shows annual mean column amounts (left) and annual and zonal 
mean mass mixing ratios (right) of cloud condensate (CO2 and H2O) for the 
same simulation. As can be seen, CO2 clouds form at high altitudes and vary 
relatively little with latitude. H2O clouds form much lower in locations that are 
dependent on the surface water sources. Because of the Martian north-south 
topographic dichotomy, H2O cloud content was much greater in the northern 
hemisphere, where temperatures were typically warmer in the low atmosphere. 
We found typical H2O cloud particle sizes of 2-10 /xm in our simulations, 
although these values were dependent on our choice of cloud condensation 
nuclei parameter [CCN] (set to 10^ kg~^ for all simulations described here; 
see Section 2 and Table 1). 



To assess the radiative impact of the H2O clouds, some simulations were per- 
formed with the H2O cloud opacity set to zero after climate equilibrium had 
been achieved. As a) global mean temperatures were low in our simulations 
even at one bar CO2 pressure and b) the planetary albedo was already sub- 
stantially modified by higher- altitude CO2 clouds, only small changes (less 
than 5 K) in the climate were observed. As we used a simplified precipitation 
scheme and assumed that the water clouds had 100% coverage in each hori- 
zonal grid cell of the model, we may have somewhat inaccurately represented 
their radiative eflPects. However, mean surface temperatures were already well 
below 273 K in our saturated water vapour simulations of Section 3.1 at the 
same pressures, so this is unlikely to have qualitatively affected our results. 
For a detailed discussion of the coverage and microphysics of CO2 clouds in 



the simulations, refer to Forget et al. [2012 



3.2.3 Transient melting 

Figure 12 (top) shows a contour plot of the annual maximum surface liquid 
H2O at each gridpoint after 40 years simulation time for the case shown in 
Figure 4 (1 bar, 25° obliquity). Given the relatively low spatial resolution 
of our model and the lack of accurate parameterisations for important sub- 
gridscale processes (slope effects, small-scale convection etc.), these results can 
only give an approximate guide to local melting under the simulated global 
climate. Nonetheless, it is clear that liquid water appears transiently in some 
amounts across the planet. The increased melting in the northern hemisphere 
can be explained by the higher temperatures there due to the same adiabatic 
effect responsible for the migration of ice to the southern highlands. 

Figure 12 (bottom) shows surface temperature against time for the four loca- 
tions shown in Figure 15. The dramatic difference in temperatures between the 
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Tharsis bulge (A) and the bottom of Hellas basin (C) is clear, along with the 
increased seasonal variation away from the equator (C and D). The transient 
melting occurring at location D (Utopia Planitia) is clear from the peak of 
surface temperatures at 273 K there between ^ 70 and 120°. Further heat- 
ing does not occur in this period because some surface H2O ice is still present. 
In Figure 13, the same results are plotted for the (f) = 45° case. There, the 
differences in insolation and absence of ice in the north mean that the great 
majority of melting events occur south of the equator. As can be seen from 
the surface temperature plot (bottom) , conditions at the equator become cold 
enough for seasonal CO2 ice deposits to form on the Tharsis bulge (location 
A; blue line). 

It is well known that Mars' obliquity evolves chaotically, and during the 
planet's history, relatively rapid changes in a range from to 70° are likely to 
have occurred many times [Laskar et al. , 2004| . With this in mind, we inves- 



tigated the effects of decreasing obliquity to zero, but keeping the end-state 
surface ice distribution of Figure 4. An obliquity of zero maximises insolation 
at the equator, which might be expected to maximise melting in the valley 
network region. However, for this case we found that while mean equatorial 
temperatures were slightly higher, the lack of seasonal variation reduced an- 
nual maximum temperatures in most locations. Hence there was a net decrease 
in the annual maximum surface liquid water (Figure 14). 

To test the sensitivity of the amount of transient melting to our assumed 
parameters, we also performed some simulations where we varied the H2O 
ice surface albedo and thermal inertia. Starting from the 1 bar, 25° obliquity 
simulation just described, we a) reduced ice albedo from 0.5 to 0.3 and b) 
increased the surface ice thermal inertia from I = 250 to 1000 J m~^ s~^/^ 
at all soil depths. In both cases all other parameters were kept constant. Figure 
16 shows the results for these two tests in terms of the annual maximum surface 
liquid water. 

Perhaps unsurprisingly, when the ice albedo was reduced melting occurred 
over a wider range of surface locations, and the total amount of melting in 
each year slightly increased. As a variety of processes may influence this pa- 



rameter, including dust transport and volcanic ash emission [e.g., Wilson and 



Head, 2007, Kerber et al., 2011 , this may have interesting implications for fu- 
ture study. However, increasing the surface ice thermal inertia essentially shut 
down transient melting entirely (Figure 16, right). As this calculation was per- 
formed assuming constant thermal inertia with depth, it clearly represents a 
lower limit on potential melting. Nonetheless, it indicates that surface seasonal 
(as opposed to basal) melting on early Mars under these conditions would al- 
most certainly be limited to small regions on the edges of ice sheets only. To 
constrain the values of these parameters better, a more detailed microphysical 
/ surface model that included the effects of sub-gridscale ice deposition, dust 
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and possibly volcanism would be required. 



3.2.4 Southern polar ice and the Dorsa Argentea Formation 



In our simulations, we find that under a dense atmosphere, thick ice sheets 
form over Mars' southern pole. Ice migrates there preferentially from the north 
because of the same adiabatic cooling effect responsible for deposition over 
Tharsis and the equatorial highlands. The Dorsa Argentea Formation is an 
extensive volatile-rich south polar deposit, with an area that may be as great 
as 2 % of the total Martian surface [Head and Pratt , 2001 . It is believed to 
have formed in the late Noachian to early Hesperian era [Plant et al., 1988 



when the atmosphere is likely to have still been thicker than it is today. Given 
that our simulations span a range of CO2 pressures, it is interesting to compare 
the results with geological maps of this region. 

As an example. Figure 17 shows surface ice in simulations at 1 bar (25 and 45° 
obliquity) and 0.2 bar (25° obliquity), with a geologic map plotted alongside 
to show the extent of the main geological features Head and Pratt, 2001 . In 



Figure 17d, yellow and purple represent the Dorsa Argentea Formation lying 
on top of the ancient heavily cratered terrain (brown) and below the current 
(Late Amazonian) polar cap (gray and white). As can be seen, there is a south 
polar H2O cap in all cases, although its latitudinal extent is greatest at 1 bar 
(Figure 17a). At current obliquity and between 0.2 and 1 bar (Figure 17a, c), 
the accumulation generally covers the area of the Dorsa Argentea Formation, 
and also extends toward 180 degrees longitude in the 0.2 bar case, a direction 
that coincides with the southernmost development of valley networks |Fassett] 
and Head, 2008b, Hynek et al., 2010 . Further simulations of ice formation 



in this region for moderately dense CO2 atmospheres would be an interesting 
topic for future research. In particular, it could be interesting to study southern 
polar ice evolution using a zoomed grid or mesoscale model, with possible 
inclusion of modifying effects due to the dust cycle. 



4 Discussion 



In contrast to the 'warm, wet' early Mars envisaged in many previous stud- 
ies, our simulations depict a cold, icy planet where even transient diurnal 
melting of water in the highlands can only occur under extremely favourable 
circumstances. As abundant quantities of liquid water clearly did fiow on early 
Mars, other processes besides those included in our model must therefore have 
been active at the time. Many cold-climate mechanisms for Martian erosion 
in the Noachian have been put forward previously, including aquifer recharge 
by hydrothermal convection [Squyres and Kasting , 1994| and flow at sub-zero 
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temperatures by acidic brines [Gaidos and Marion, 2003, Fairen, 2010 . In 
the following subsections, we focus on three particularly debated processes: 
heating by impacts, volcanism, and basal melting of ice sheets. 



4.1 Meteorite impacts 



As noted in Section 1, impacts have already been proposed by [Segura et al. 



2002, 2008 



as the primary cause of the valley networks via the formation of 
transient steam atmospheres. However, one of the key Segura et al. arguments, 
namely that the rainfall from these atmospheres post-impact is sufficient to 
reproduce the observed erosion, has been strongly criticized by other authors 
[Ramirez et al. , 2012 . Alternatively, one could imagine that the numerous im- 
pacts that occurred during the Noachian (including smaller ones) would have 
induced local melting of surface ice if glaciers were present. Such a mechanism 
has been suggested to explain fluvial landforms on fresh Martian impact eject a 
[Mangold , 2012 . Under the higher atmospheric pressure and greater impact or 
flux of the Noachian, it could have caused much more extensive fluvial erosion. 
However, for impact-induced melting to be a viable explanation for the valley 
networks, some mechanism to preferentially transport water to the southern 
highlands must be invoked. 



In our simulations, under the denser CO2 atmosphere that would be ex- 
pected due to Tharsis rise formation, ice was deposited and stabilized in the 
Noachian highland regions due to the adiabatic cooling effect. Under these cir- 
cumstances, heating due to impacts could cause extensive melting and hence 
erosion in exactly the regions where the majority of valley networks are ob- 
served Fassett and Head, 2008b, Hynek et al., 2010 . Transient melting would 



transport water to lower lying regions, but once temperatures again dropped 
below the freezing point of water, the slow transport of ice to the highlands 
via sublimation and snowfall would recommence. As in our simulations that 
started with ice at low altitudes, the climate system would then return to a 
state of equilibrium on much longer timescales. Figure 18 shows a schematic 
of this process. 



Although our results provide a possible solution to one of the problems as- 
sociated with impact-dominated erosion, the impact hypothesis remains con- 
troversial in light of some other observations. For example, young, very large 
impact craters, such as the Amazonian-era Lyot, have few visible effects of re- 
gional or global-scale melting. This is often regarded as further evidence that 
impacts were minor players in valley network erosion in the Noachian. Fur- 
ther study of hydrology and climate under the extreme conditions expected 
post-impact are therefore needed to assess the plausibility of this scenario in 
detail. 
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4^2 Volcanism 



Another potential driver of climate in the Noachian that we have neglected 
in these simulations is volcanism. As well as causing substantial emissions 
of CO2, late-Noachian volcanic activity may have influenced the climate via 
emission of sulphur gases (SO2 and H2S) and pyroclasts (dust/ash particles) 
[Halevy et al. , 2007, Kerber et al., 2012| . One-dimensional radiative-convective 
studies have estimated the sulphur gases to have a potential warming effect 
of up to 27-70 K |Postawko and Kuhn, 1986, Tian et al., 2010 , while dust is 
most likely to cause a small (2-10 K) amount of warming, depending on the 



particle size and vertical distributions | Forget et al., 2012| . 



As mentioned in the introduction, Tian et al. [2010 argued that the net ef- 



fect of volcanism should be to cool the early climate, due to the rapid aerosol 
formation and hence anti-greenhouse effect that could occur after sulphur 
gas emission. However, their conclusions were based on the results of one- 
dimensional climate simulations performed without the effects of CO2 clouds 
included. As described in Forget et al. 2012 , carbon dioxide clouds have a 



major impact on the atmospheric radiative budget by raising planetary albedo 
and increasing downward IR scattering. It is therefore likely that aerosol ra- 
diative effects would be substantially different in a model where clouds were 
also included. Furthermore, sulphate aerosol particles in the C02-condensing 
region of an early Martian atmosphere would be potential condensation nuclei 
for the formation of the CO2 ice cloud crystals, which would further influ- 
ence the distribution of both aerosol and cloud particles. An extension of the 
work presented here that included the chemistry and radiative effects of sul- 
phur compounds would therefore be an extremely interesting avenue of future 
research. 



4-3 Basal melting 



It has been proposed that at least some Noachian fluvial features could have 



formed as a result of the basal melting of glaciers or thick snow deposits [Carr 



1983, Carr and Head, 2003]. Geological evidence suggests that throughout 



much of the Amazonian, the mean annual surface temperatures of Mars have 
been so cold that basal melting does not occur in ice sheets or glaciers [e.g., 
Fastook et aTj |2008[ [Head et al.[ |2010[ [Fastook et alj |2011| . However, the doc- 



umented evidence for extensive and well-developed eskers in the Dorsa Argen- 
tea Formation indicates that basal melting and wet-based glaciation occurred 
at the South Pole near the Noachian-Hesperian boundary. Recent glacial ac- 
cumulation and ice-flow modeling Fastook et al., 20121 has shown that to 



produce signiflcant basal melting for typical Noachian-Hesperian geothermal 



18 



heat fluxes (45-65 mW m~^), mean annual south polar atmospheric tempera- 
tures of -50 to -75 °C are required (approximately the same range as we flnd in 
our 1 bar simulations). If geologically based approaches to constraining south 
polar Noachian/Hesperian temperatures such as this prove to be robust, this 
provides further evidence in favour of a cold and relatively dry early Martian 
climate. 



4^4 Conclusions 



Our results have shown that early Mars is unlikely to have been warm and wet 
if its atmosphere was composed of CO2 and H2O only, even when the effects 
of CO2 clouds are taken into account. It is possible that other greenhouse 
gases or aerosols due to e.g. Tharsis volcanic emissions also played a role in 
the climate. However, given current constraints on the maximum amount of 



CO2 in the primitive atmosphere (e.g., Grott et al. 2011 ; see also Forget 



et al. 2012j for a detailed discussion), it seems improbable that there was any 



combination of effects powerful enough to produce a steady-state warm, wet 
Mars at any time from the mid-Noachian onwards. Despite this, in climates 
where global mean temperatures are below zero there are still effects that 
can contribute to transient melting and hence erosion. Simulating the water 
cycle in 3D has shown that even when the planet is in a state of thermal 
equilibrium, with mean temperatures as low as ^230 K, seasonal and diurnal 
warming can cause some localised melting of H2O, although the amount that 
can be achieved depends strongly on the assumed surface albedo and thermal 
inertia. 

As the era of Tharsis formation / volcanism (and hence of increased CO2 at- 
mospheric pressure) that we have modelled here occurred at a time of elevated 
meteorite bombardment, impacts will also have repeatedly caused melting of 
stable ice deposits in the Noachian highlands. Even if post-impact rainfall 
was relatively short-term, the fact that ice always returned to higher terrain 
due to the adiabatic cooling effect means that each impact could have cre- 
ated signiflcant amounts of meltwater in the valley network regions. Volcanic 
activity, particularly that associated with formation of the Tharsis bulge and 



Hesperian ridged plains Head et al., 2002 , could also have caused discrete 



warming events in theory, although it remains to be demonstrated how the 
anti-greenhouse effect due to sulphate aerosol formation could be overcome in 
this scenario. 

Although we believe these simulations have probably captured the main global 
features of the steady-state Martian climate in the late Noachian, there are 
clearly a number of other interesting effects that could be investigated in 
future. We have not included dust or volcanic emissions in our simulations. 
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Both of these are hkely to have had an important effect during the Noachian, 
and adding them in future would aUow a more complete assessment of the 
nature of the early water cycle. Modelling the effects of impacts directly in 
3D could also lead to interesting insights, although robust parameterization of 
many physical processes across wide temperature ranges and short timescales 
would be necessary to do this accurately. A more immediate application could 
be to couple the climate model described here with a more detailed subsurface 
/ hydrological scheme. Such an approach could be particularly revealing in 
localised studies of cases where geothermal processes or residual heating due to 
impacts may have been important [e.g., Abramov and Kring, 2005| . In future, 
integration of 3D climate models with specific representations of warming 
processes (both local and global) should allow a detailed assessment of whether 
the observed Noachian fluvial features can be reconciled with the cold, mainly 
frozen planet simulated here. 
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Fig. 1. Comparison of present-day Martian topography (top) with the reduced Thar- 
sis bulge topography (bottom) described by equation (4) and used in some simula- 
tions. 
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Fig. 2. Effects of atmospheric CO2 and H2O on global temperature. Error bars show 
mean and maximum / minimum surface temperature vs. pressure (sampled over one 
orbit and across the surface) for dry CO2 atmospheres (red), and simulations with 
100% relative humidity (blue) but no H2O clouds. Dashed and dotted black lines 
show the condensation curve of CO2 and the melting point of H2O, respectively. 
For this plot simulations were performed at 0.2, 0.5, 1 and 2 bar; the dry and wet 
data are slightly separated for clarity only. 
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Fig. 3. Surface temperature mean and variations in the 100% relative humidity 
simulations. Top, middle and bottom plots show the annual mean, the diurnally 
averaged annual maximum, and the absolute annual maximum, respectively. Left, 
right and center columns are for 0.008, 0.2 and 1 bar, respectively. Black contours 
show the topography used in the simulation. Strong correlation of annual mean 
surface temperature with altitude (as occurs on e.g. present-day Earth and Venus) 
is apparent in the 0.2 and 1 bar cases. 
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Fig. 4. Evolution of H2O ice in simulations using the iteration scheme described in 
Section 2. From top to bottom, plots show surface ice density in kg m~^ at the start 
of the simulation and annual mean after 4, 20 and 40 Mars years. Ice iteration was 
performed every two years, with a 100-year timestep used for the first 5 iterations 
and 10-year timesteps used thereafter. Simulations were performed at 1 bar mean 
surface pressure with obliquity 25°. Left and right columns show cases with initial 
ice deposits in low-lying regions and at the poles, respectively. 
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Fig. 5. Zonally and yearly averaged H2O surface ice density as a function of latitude 
and time for the same simulations as in Figure 4. The dotted vertical line indicates 
the transition from 100-year to 10-year timesteps in the ice evolution algorithm. 
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Fig. 6. Same as Figure 4 (right column) except with modified topography (left) and 
with obliquity increased to 45° (right). 
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Fig. 7. Same as Figure 4 (left column) except for surface pressure of 0.04 bar (left) 
and 0.008 bar (right). 
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Fig. 8. Annual mean horizontal wind at 500 hPa (left) and annual and zonal mean 
vertical velocity (right), for the 1-bar 25° obliquity case after 40 simulation years. 
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Fig. 9. Same as Figure 8, but for 45° obliquity. 
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Fig. 10. Total precipitation (snowfall) in each season for Lg = 0° to 90°, 90° to 180°, 
180° to 270° and 270° to 360° (in descending order), for the 1-bar 25° obliquity case 
after 40 simulation years. 
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Fig. 11. Annual mean column amounts (left) and annual and zonal mean mass 
mixing ratios (right) of CO2 (red) and H2O (blue) cloud condensate for the 1-bar 
25° obliquity case after 40 simulation years. 
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Fig. 12. (top) Maximum surface liquid water in one year after 40 years simulation 
time, for the 1-bar simulation with obliquity 25° shown in Figure 4. (bottom) Surface 
temperature vs. time for the same simulation, for the four locations displayed in 
Figure 15. 
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Fig. 13. Same as Figure 12, except for obliquity 45°. 
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Fig. 14. Same as Figure 12, except with obliquity set to 0° after 40 years simulation 
time. 



80' 



7 J \ ' 








D _ - 

\ 


r 


R 


, r 







180°W120°W 60°W 0° 60°E 120^E 180°W 

Fig. 15. Locations for the four temperature time series displayed in Figs. 12-14. 
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Fig. 16. Maximum surface liquid water in one year after 40 years simulation time, 
for the 1-bar simulation with obliquity 25° shown in Figure 4, after (left) change of 
ice surface albedo from 0.5 to 0.3 and (right) increase of ice surface thermal inertia 
from 250 to 1000 J m-^ g-^s k-^. 
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Fig. 17. Plots of the annual mean water ice coverage on the southern pole after 40 
years for mean surface pressure 1 bar and obliquity a) 25°, b) 45° and mean surface 
pressure 0.2 bar and obliquity 25°. d) For comparison, a plot of the Dorsa Argentea 
Formation fromlHead and Prattl 120011 is also shown. 
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Fig. 18. Schematic of the effect of periodic melting events under a moderately dense 
CO2 atmosphere, a) In a steady state, ice deposits are concentrated in the colder 
highland regions of the planet, b) Impacts or volcanism cause transient ice melting 
and flow to lower lying regions on short timescales. c) On much longer timescales, 
ice is once again transported to highland regions via sublimation and light snowfall. 



37 



